Setup

# Load packages here
library(tidyverse)
library(factoextra)
library(gridExtra)     
library(sf)
library(factoextra)
library(ggrepel)

Data

Load the file countries.Rdata.

load("countries.Rdata")

It contains two objects: df is a data frame of 138 countries with information on 10 variables (see short descriptions below). sf is a simple features object that contains the geometries of the countries’ borders.

life_expectancy: average number of years a newborn child is expected to live
HDI: index that ranks countries by level of human development in terms of health, education, and living standard
income_person: GDP per capita
gini_coefficient: income inequality - a higher number means more inequality.
water: percentage of people using at least basic water services.
sanitation: percentage of people using at least basic sanitation services
calories: measures the energy content of the food.
freedom: index of political rights and civil liberties, on a range from 1 (most free) to 7 (least free)
democracy: index of quality of democracies between 0 and 100
corruption: score of perceptions of corruption by Transparency International. From 0 (highly corrupt) to 100 (very clean).
broadband: fixed subscriptions to high-speed access to the public Internet
internet_users: internet users in percentge of population
covid_confirmed: number of confirmed covid cases until 2021-06-05 per 1000 inhabitants

Exercise 1: Clustering

Carry out a hierarchical clustering analysis.

Exercise 1.1

Check whether it is necessary to preprocess the data. Justify your choice, and (if applicable) explain what the preprocessing steps are.

Antwort

Ja, es ist notwendig, die Daten vorher zu verarbeiten. Zum Einen haben die Daten unterschiedliche Einheiten. Wir haben es mit der Lebenserwartung in Jahren zu tun und dem prozentualen Anteil der Bevölkerung, die einen Zugang zu Wasser hat. Zudem die Menge an Kalorien, die eine Person am Tag zu sich nimmt, etc.. Dazu kommt, dass die einzelnen Spalten unterschiedliche Datenweiten haben. Der Vergleich zwischen den Spalten ist somit nicht möglich, wir würden sonst Äpfel mit Birnen vergleichen. Die Spalte “iso_alpha” ist außerdem nicht numerisch.

Die Daten müssen also zunächst standardisiert werden. Die Standardisierung geschieht, indem von den Daten die Mittelwerte der jeweiligen Spalten abgezogen werden, sodass der Mittelwert 0 ist. Danach müssen die Daten durch die Standardabweichungen geteilt werden.

df %>%
  head()

df %>%
  select(-iso_alpha) %>%
  apply(2, mean)

##  life_expectancy              HDI    income_person gini_coefficient 
##     7.269118e+01     7.044706e-01     1.763697e+04     3.870809e+01 
##            water       sanitation         calories          freedom 
##     8.659926e+01     7.319941e+01     2.876471e+03     3.246324e+00 
##        democracy       corruption        broadband   internet_users 
##     5.915809e+01     4.417647e+01     1.256747e+01     4.849096e+01 
##  covid_confirmed 
##     3.603676e+01

df %>%
  select(-iso_alpha) %>%
  apply(2, sd)

##  life_expectancy              HDI    income_person gini_coefficient 
##     7.730195e+00     1.591161e-01     1.737466e+04     8.169873e+00 
##            water       sanitation         calories          freedom 
##     1.656978e+01     2.939221e+01     4.637616e+02     1.881682e+00 
##        democracy       corruption        broadband   internet_users 
##     2.060261e+01     1.968714e+01     1.326526e+01     2.777444e+01 
##  covid_confirmed 
##     3.632373e+01

Exercise 1.2

Name three possible distance measures for hierarchical clustering. Select an appropriate distance measure for the given type of data and justify your choice. Visualize a distance matrix for it.

Antwort

Drei mögliche Distanzmaße für das hierarchische Clustering sind die euklidische, die Kosinus und die Pearson Distanz. Die für diesen Datensatz sinnvollste Distanz ist die euklidische Distanz, da wir die Staaten beispielsweise in verschiedene Entwicklungsstufen einteilen wollen. Dafür müssen wir die konkreten Abstände zwischen den Werten der Variablen im DataFrame wissen. Bei der Pearson Distanz könnten wir eher erkennen, ob sich Länder in bestimmten Bereichen ähnlich sind, also zum Beispiel, dass Argentinien in jeder Kategorie des Datensatzes einen doppelt so hohen Wert hat, wie Portugal (Nur ein Beispiel). Die Pearson Distanz würde die beiden Länder dann mit einer geringen Distanz einstufen.

# The distance matrix will be large, adjust the fig.height such that is can be acceptably recognized

d_eucl <- df %>%
  select(-(iso_alpha)) %>%
  get_dist(method = "euclidean", stand = TRUE)

fviz_dist(d_eucl) +
  geom_text(aes(label = round(value, 2)), size = 2)

Exercise 1.3

Run agglomerative clustering with 4 different linkage methods and plot the dendrogramms. Assess the usefulness of these linkage methods, and choose your preferred linkage method. Justify your choice.

Anwort

Die Linkage Methoden unterscheiden sich darin, auf welche Art und Weise zusammengehörende Cluster bestimmt werden. 1. Complete Linkage: Die größte Distanz zwischen zwei Datenpunkt verschiedener Cluster wird berechnet. 2. Single Linkage: Die kleinste Distanz zwischen zwei Datenpunkt verschiedener Cluster wird berechnet. 3. Average Linkage: Die durchschnittliche Distanz zwischen zwei Datenpunkt verschiedener Cluster wird berechnet. 4. Ward’s Minimum variance criterion: Hierbei werden in jedem Schritt die beiden Cluster miteinander verbunden, die zu dem geringsten Anstieg in der Varianz führen.

Die am wenigsten sinnstiftende Methode ist hier die single Linkage Methode. Man sieht, dass mit jedem weiteren Schritt innerhalb des Dendogramms nur ein weiteres Land geclustert wird. Diese Methode dient eher der Abgrenzung eines einzelnen Landes, als der gröberen Gruppierung mehrerer Staaten.

Die Average Linkage Methode hat bei vier Clustern interessanterweise 3 größere Cluster hervorgebracht und ein kleines Cluster, das drei der 7 Staaten der Arabischen Halbinsel umfasst.

Die Complete Linkage und Wards Methode haben beide ungefähr gleich große Ergebnismengen hervorgebracht. Die Complete Linkage Methode legt aber einen größeren Fokus darauf, stark unterschiedliche Datenpunkte voneinander zu unterscheiden, wodurch wahrscheinlich eine bessere Gruppierung erreicht wird, als nur darauf zu achten, wie stark sich die Datenpunkte in einem einzelnen Cluster unterscheiden.

Daher ist meine bevorzugte Clustering-Methode die Complete Linkage Methode, da diese mir sowohl in der Theorie, als auch in den Ergebnissen augenscheinlich logischere Gruppierungen ergibt.

# The dendrogram will be large, adjust the fig.height such that is can be acceptably recognized
result_complete <- hcut(
    x = select(df, -iso_alpha), 
    k = 4,
    hc_func = "agnes",          
    hc_metric = 'euclidean',  #with euclidean a histogram
    hc_method = 'complete',
    stand = TRUE)
  fviz_dend(result_complete, horiz = TRUE, cex = 0.8, labels_track_height =  3,
            main = glue::glue("Method: complete, Distance: euclidean"))

# The dendrogram will be large, adjust the fig.height such that is can be acceptably recognized
result <- hcut(
    x = select(df, -iso_alpha), 
    k = 4,
    hc_func = "agnes",          
    hc_metric = 'euclidean',  #with euclidean a histogram
    hc_method = 'single',
    stand = TRUE)
  fviz_dend(result, horiz = TRUE, cex = 0.8, labels_track_height =  3,
            main = glue::glue("Method: single, Distance: euclidean"))

# The dendrogram will be large, adjust the fig.height such that is can be acceptably recognized
result <- hcut(
    x = select(df, -iso_alpha), 
    k = 4,
    hc_func = "agnes",          
    hc_metric = 'euclidean',  #with euclidean a histogram
    hc_method = 'average',
    stand = TRUE)
  fviz_dend(result, horiz = TRUE, cex = 0.8, labels_track_height =  3,
            main = glue::glue("Method: average, Distance: euclidean"))

# The dendrogram will be large, adjust the fig.height such that is can be acceptably recognized
result <- hcut(
    x = select(df, -iso_alpha), 
    k = 4,
    hc_func = "agnes",          
    hc_metric = 'euclidean',  #with euclidean a histogram
    hc_method = 'ward.D2',
    stand = TRUE)
  fviz_dend(result, horiz = TRUE, cex = 0.8, labels_track_height =  3,
            main = glue::glue("Method: Ward.D2, Distance: euclidean"))

Exercise 1.4

Analyse how many clusters are suitable (in the sense of being informative, interesting, distinguishable), based on:

visual inspection of the dendrogramms
summary statistics (e.g. average characterists of each cluster)
A map of the world, where each country is colored according to its cluster affiliation

Justify, how many clusters you would choose. Explain what characterises your clusters. And also mention if there are cluster affiliations which are surprising to your eyes.

Antwort

Meine bevorzugte Wahl bleibt bei den 4 Clustern. Anhand des Dendogramms ist zu erkennen, dass die nächste sinnvolle Clusteranzahl bei einem deutlich höheren Wert als 4 liegen würde. Geht man in dem Dendogramm nämlich tiefer, spalten sich annähernd gleichzeitig alle 4 Cluster ein weiteres Mal,sodass dann 8 Cluster entstehen würden. Dies erscheint mir aber nicht sinnvoll, so viele Cluster zu erzeugen.

Bei der Berechnung der Mittelwerte innerhalb eines Clusters zeigt sich zudem, dass sich die Unterschiede zwischen den Clustern in den Kategorien Lebenserwartung, HDI, Einkommen und Internet-Users am offensichtlichsten bei 4 Clustern zeigen. Wenn man sich die Cluster-Mittelwerte ansieht, und die Cluster nach der Lebenserwartung sortiert, erkennt man, dass diese Sortierung auch bei fast allen anderen Variablen zu erkennen ist. Nur beim gini_coefficient, democracy und freedom sind die Cluster an zweiter und dritter Stelle vertauscht.

Auffällig ist bei der Einfärbung der Weltkarte die offensichtlich gleiche Einfärbung von Europa, USA, Kanada und Australien.

result_complete <- hcut(
    x = select(df, -iso_alpha), 
    k = 4,
    hc_func = "agnes",          
    hc_metric = 'euclidean',  #with euclidean a histogram
    hc_method = 'complete',
    stand = TRUE)
  fviz_dend(result_complete, horiz = TRUE, cex = 0.8, labels_track_height =  3,
            main = glue::glue("Method: complete, Distance: euclidean"))

df %>%
  select(-iso_alpha)%>%
  mutate(cluster = result_complete$cluster) %>%
  group_by(cluster)%>%
  summarize(across(everything(), mean)) %>%
  arrange(by = desc(life_expectancy))

df %>%
  left_join(sf)%>%
  sf::st_as_sf()%>%
  mutate(cluster = result_complete$cluster) %>%
  ggplot()+
  geom_sf(aes(fill = factor(cluster)))+
  labs(title = "Cluster laut der complete Linkage Methode mit 4 Clustern")

### Nur als Spielerei noch den HDI als Kategorie, wie die vereinten Nationen ihn kategorisieren.
df %>%
  left_join(sf)%>%
  sf::st_as_sf()%>%
  mutate(cluster = case_when(
    HDI > 0.8 ~ 1,
    HDI > 0.7 ~ 2, 
    HDI > 0.55 ~3,
    TRUE ~ 4
  )) %>%
  ggplot()+
  geom_sf(aes(fill = factor(cluster)))+
  labs(title = "HDI-Kategorien laut der UN")

Exercise 2: PCA

Exercise 2.1

If necessary preprocess the data. Then carry out a principal component analysis and show a biplot.

Explain which patterns can be recogized from the biplot, using the following countries and variables as examples (Please do not comment on any combination of country and variable, only on selected informative patterns):

countries: Norway, Chad, Afghanistan, South Africa, Peru, Botswana
variables: democracy, gini_coefficient, life_expectancy, covid_confirmed

Antwort

Zunächst kann man erkennen, dass die erste PC mit fast allen Variablen aus dem Datensatz positiv korreliert. Einzig der gini_coefficient und freedom sind negativ korreliert und stellen neben der covid_confirmed Variable die einzigen beiden Variablen dar, die eher negativ behaftet sind, wenn sie einen hohen Wert aufweisen (ein hoher gini_coefficient-Wert steht für hohe Ungleichheit und ein hoher freedom Wert steht für wenig Freiheit). Im Gegensatz dazu sind alle anderen Variablen mit hohen Werten eher positiv behaftet (ein hoher democracy Wert steht für bessere Demokratie, eine hohe life_expectancy steht für ein langes Leben). Insofern können wir davon ausgehen, dass Länder, die einen hohen positiven Score auf der ersten Principal Komponent haben, einen hohen Lebensstandard ermöglichen. Länder, die dagegen einen negativen Score haben, sind eher lebensunwürdig und bieten weniger essentielle Lebensgrundlagen.

Interessanterweise korrelieren die democracy und die freedom Variable fast exakt negativ miteinander. Die gini_coefficient Variable ist dagegen gänzlich unkorreliert mit der democracy und der freedom Variable, was dafür spricht, dass die Demokratie eines Landes nichts mit der Ungleichheit zu tun hat. Außerdem fällt auf, dass die meisten Variablen ungefähr gleich stark auf die beiden Principal Components “laden”. Die stärksten Variablen sind aber die democracy und die freedom Variable, während die covid_confirmed die schwächste Variable ist.

Einzelne Variablen: - democracy: Laut dem Democracy loading, zählt Norwegen zu den am stärksten demokratisierten Ländern, während Chad auf einen der letzten Plätze fällt. - gini_coefficient: Botswana und South Africa zählen zu den Ländern mit dem am stärksten ungleich verteilten Vermögen, während in Norwegen das Vermögen gleicher verteilt zu sein scheint. - life_expectancy: Die Lebenserwartung korreliert stark mit dem prozentualen Zugang zu Wasser, der Kalorienaufnahme und dem Zugang zu sanitären Einrichtungen. Länder wie Belarus und Kuwait haben hier einen ähnlich hohen Score wie Japan und Australien. In Ländern wie Chad und Afghanistan ist die Lebenserwartung sehr niedrig. - covid_confirmed: Die covid_confirmed-Variable korreliert sehr stark mit der ersten PC und damit auch mit dem income_person. Länder, die also einen höheren Lebensstandard haben, scheinen somit auch stärker von Covid-19 betroffen zu sein.

Peru scheint auf allen einzelnen Variablen und den PCs einen niedrigen score zu besitzen. Es ist also ein stark mittelmäßiges Land, das keine besondere Eigenschaft bei den einzelnen Variablen hat.

df_pca <- df %>% select(-iso_alpha) %>% scale()
pca <- prcomp(df_pca, center = FALSE, scale.=FALSE)
fviz_pca_biplot(pca, 
                repel = TRUE, 
                geom = "text" 
                )

Exercise 2.2

Visualize the cummulative sum of the percentage of variance explained (PVE) as a function of the number of pincipal components. What is the fraction of variance explained by the first two principle components? And based on this answer, how do you judge the explanatory power of the biplot: does the biplot give us a meaningful 2-dimensional representation of the patterns in the data?

Antwort

Die ersten beiden Principal Components erklären etwa 77.6% der gesamten Varianz in dem Datensatz, wovon die erste PC bereits 67.6% erklärt. Der Biplot gibt uns somit bereits eine sehr bedeutsame Repräsentation des Datensatzes. Diese hohe PVE-Wert kommt wahrscheinlich auch daher, dass sehr viele der Variablen in dem Datensatz bereits miteinander korrelieren und somit viele Informationen durch weitere Variablen nur redundant sind. Die Lebenserwartung und der Zugang zum Wasser beispielsweise sind beinahe 100% korreliert, womit eine der beiden Variablen auch weg gelassen werden könnte und der PVE-Wert nur sehr wenig geringer würde.

scores <- pca$x


var <- cov(scores) %>% diag()

stats <- tibble(pc = factor(1:13),        
           variance = var,                    
           pve = var / sum(var),               
           cum_pve = cumsum(pve))               
stats %>%
  ggplot() + 
  geom_col(aes(pc, cum_pve, fill = "Cummulative PVE")) +
  geom_col(aes(pc, pve, fill = "PVE")) +
  labs(x = "Principle Components", 
       y= "Percentage of Variance Explained",
       fill = "Type", 
       title = "Erklärte Varianz unter Verwendung entsprechender PC")

stats$cum_pve[2]

##       PC2 
## 0.7759666

Exercise 2.3

Keep only the first two principle components, and discard the other ones. Then approximately re-create the original data set (unscaled and uncentered) using only these two principal components. Briefly state how well the approximation works, based on rough inspection of the true and the approximated data.

Ansatz

Wir können den ursprünglichen Datensatz wiederherstellen, indem wir die Scores aus den ersten beiden Spalten der Scores-Matrix mit der Transponierten Weight-Matrix(den ersten beiden Spalten dieser) multiplizieren. Danach müssen wir die Daten nur wieder mit der Standardabweichung der originalen Daten multiplizieren und den Mittelwert addieren.

Die Approximation funktioniert relativ gut, es sind bei einigen Variablen aber trotzdem starke Abweichungen zu erkennen. Zudem gibt es in dem approximierten Datensatz auch Variablen in einem Wertebereich, der gar nicht exitiert, zum Beispiel hat Afghanistan ein negatives income_person oder Luxemburg einen Wasserzugang von über 112%.

Die Histogramme zeigen auch Übereinstimmungen, aber es sind auch deutliche Unterschiede zu sehen. Die beiden ersten Principal Components erklären eben nur knapp 78% der Varianz, dementsprechend stimmt der zurückgewonnene Datensatz auch nur zu 78% dem Originaldatensatz.

scores <- pca$x[,1:2 ]
weight <- pca$rotation
sd <- pca$scale


recovered_data <- t(t(pca$x[,1:2 ] %*% t(pca$rotation[,1:2 ])) * pca$scale + pca$center) %>%
  as.data.frame()

recovered_data %>%
  summary()

##  life_expectancy      HDI         income_person    gini_coefficient
##  Min.   :57.89   Min.   :0.3805   Min.   :-12959   Min.   :28.95   
##  1st Qu.:66.37   1st Qu.:0.5721   1st Qu.:  5412   1st Qu.:34.64   
##  Median :73.89   Median :0.7076   Median : 16871   Median :38.01   
##  Mean   :72.69   Mean   :0.7045   Mean   : 17637   Mean   :38.71   
##  3rd Qu.:78.66   3rd Qu.:0.8374   3rd Qu.: 30044   3rd Qu.:42.78   
##  Max.   :85.40   Max.   :0.9980   Max.   : 46308   Max.   :50.62   
##      water          sanitation        calories       freedom         
##  Min.   : 56.36   Min.   : 17.31   Min.   :2032   Min.   :-0.003402  
##  1st Qu.: 73.61   1st Qu.: 50.34   1st Qu.:2532   1st Qu.: 1.807861  
##  Median : 89.08   Median : 77.81   Median :2947   Median : 3.228362  
##  Mean   : 86.60   Mean   : 73.20   Mean   :2876   Mean   : 3.246324  
##  3rd Qu.: 98.79   3rd Qu.: 96.83   3rd Qu.:3232   3rd Qu.: 4.657657  
##  Max.   :112.51   Max.   :118.82   Max.   :3567   Max.   : 6.648098  
##    democracy       corruption      broadband       internet_users 
##  Min.   :21.71   Min.   :11.72   Min.   :-12.257   Min.   :-7.09  
##  1st Qu.:42.79   1st Qu.:29.35   1st Qu.:  2.053   1st Qu.:26.04  
##  Median :59.36   Median :42.89   Median : 11.762   Median :49.70  
##  Mean   :59.16   Mean   :44.18   Mean   : 12.567   Mean   :48.49  
##  3rd Qu.:75.66   3rd Qu.:57.55   3rd Qu.: 22.603   3rd Qu.:71.68  
##  Max.   :96.13   Max.   :79.52   Max.   : 36.434   Max.   :98.53  
##  covid_confirmed 
##  Min.   :-19.63  
##  1st Qu.: 14.10  
##  Median : 34.96  
##  Mean   : 36.04  
##  3rd Qu.: 58.57  
##  Max.   : 87.80

Histogramm der beiden Datensätze

columns <- colnames(recovered_data)

df_orig <- df %>%
  mutate("dataframe" = "original") %>%
  select(-iso_alpha) %>%
  rownames_to_column()

combined_df <- 
  recovered_data %>%
  rownames_to_column()  %>%
  mutate("dataframe" = "recovered") %>%
  rbind(.data, df_orig) %>%
  pivot_longer(cols = columns, names_to = "variables", values_to = "value")

ggplot(combined_df, aes(value, fill = dataframe)) + 
  geom_density(alpha = 0.5, position = "identity") +
  facet_wrap(variables ~ ., scales = "free", ncol = 2)

Exercise 2.4

Calculate the loadings matrix, i.e. the correlations between the original (but possibly preprocessed) data and the principal components. And state which of the variables are not well represented by the first two principal components.

Antwort

Die am wenigsten stark von den ersten beiden PCs repräsentierten Variablen sind covid_confirmed und income_person, die jeweils nur einen loading score von 0.24 und 0.28 besitzen.

Die am stärksten repräsentierten Variablen sind democracy und freedom mit einem score von 0.58 und 0.6.

weight <- pca$rotation
data.frame(weight[, 1:2]) %>%
  mutate(eucl_length = sqrt(PC1 * PC1 + PC2 * PC2)) %>%
  arrange(eucl_length)

Then prove the following properties:

The Euclidean norm of the row vectors of the loadings matrix are equal to 1 (it is sufficient to show this for 1 arbitrarily chosen row vector)

sum_ = 0
for(v in weight[1,]){
  sum_ = sum_ + v^2
}
print(sqrt(sum_))

## [1] 1

The principle components (i.e. the columns of the scores matrix) are uncorrelated with each other

scores <- pca$x
cor(scores) %>%
  round(2) %>%
  data.frame()

The pairwise dot products of the principle components is equal to 0. What is the geometric interpretation of a dot product equal to 0? (it is sufficient to show this for 1 arbitrarily chosen dot product)

Antwort

Das Ergebnis eines Skalarproduktes zweier Vektoren gleich 0 bedeutet, dass diese beiden Vektoren aufeinander orthogonal stehen, sie sind als nicht korreliert. Im Endeffekt bedeutet das dasselbe wie die Korrelationsmatrix von der Aufgabe darüber. Das Skalarprodukt eines Vektors mit sich selbst bedeutet hingegen, dass beide Vektoren in die gleiche Richtung zeigen, wenn es sich um normierte Vektoren handelt. Die PCs sind normiert, dadurch kann dieser Schluss gezogen werden.

Geometrisch interpretiert handelt es sich also um einen Vektor, der um 90° gegenüber dem anderen Vektor gedreht ist.

weight[,2] %*% weight[,3]    # ergibt 0, also orthogonal

##              [,1]
## [1,] -2.02529e-16

weight[,2] %*% weight[,2]    # ergibt 1, also gleiche Richtung

##      [,1]
## [1,]    1

Portfolio Project 2: Clustering and PCA

Setup

Data

Exercise 1: Clustering

Exercise 1.1

Antwort

Exercise 1.2

Antwort

Exercise 1.3

Anwort

Exercise 1.4

Antwort

Exercise 2: PCA

Exercise 2.1

Antwort

Exercise 2.2

Antwort

Exercise 2.3

Ansatz

Histogramm der beiden Datensätze

Exercise 2.4

Antwort

Antwort